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d ■ Abstract. We have performed a systematic study of the effects of field strength and 

quenched disorder on the driven dynamics of square artificial spin ice. We construct a 
^ ^ _ network representation of the configurational phase space, where nodes represent the 

^ ■ microscopic configurations and a directed link between node i and node j means that 

O , the field may induce a transition between the corresponding configurations. In this 

way, we are able to quantitatively describe how the field and the disorder affect the 
connectedness of states and the reversibility of dynamics. In particular, we have shown 
(N- that for optimal field strengths, a substantial fraction of all states can be accessed using 

external driving fields, and this fraction is increased by disorder. We discuss how this 
\^ ■ relates to control and potential information storage applications for artificial spin ices. 
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A network model for field and quenched disorder effects in artificial spin ice. 2 
1. Introduction 

In an idealised model, the nanoscale magnetic islands of artificial spin ice P-[3] are 
identical Ising macrospins. However, in reality, unavoidable small variations during 
the island fabrication process lead to a distribution of island properties. In fact, a 
comparison of experimental and simulation results for square artificial spin ice shows 
that the effective distribution of island magnetisation switching barriers has a width 
on the same scale as nearest-neighbour interactions Quenched disorder provides 
pinning and nucleation sites and modifies hysteresis [SHTO]. Accordingly, a complete 
understanding of artificial spin ice requires not only an understanding of the behaviour 
of frustrated coupled Ising spins subject to an external magnetic field, but also an 
understanding of the role of quenched disorder in those dynamics. 

One approach to studying square artificial spin ice dynamics that has received 
much attention treats the vertices of the array as objects [IlSlEllTlHST] . For example, the 
populations of different vertex types (distinguished by their energy) provide a measure of 
the level of ordering and can be analysed in terms of an effective temperature associated 
with ac demagnetisation; magnetisation reversal of an array can be characterised by 
the motion and interactions of 'monopole' vertices; the effects of quenched disorder that 
affects interactions can be described in terms of variations in vertex energies; and the 
evolution under a rotating applied field can be modelled in terms of population dynamics 
of vertex types, allowing analytical expressions for the system's evolution to be written 
down and solved. 

We have recently demonstrated [22] the value of an alternative approach, in which 
the fundamental objects of interest are not vertices but whole-array spin configurations. 
In this approach, the set of all Ising spin configurations forms a discrete phase space, 
and the action of an applied field is to 'transport' the system from one point in 
phase space to another, via one or more spin flips. This picture is essentially a 
mapping of dynamics onto a directed network, in which configurational states are 
nodes and a directed hnk exists from node i to node /, i.e., z — )• /, if an applied 
field can drive the system from configuration i to configuration /. In other words, the 
network describes which barriers to flipping spins may be overcome by an external field. 
Related approaches involving network analysis have been used previously in the study of 
dynamical maps [23112^ . geometrically frustrated systems [251427] . the random field Ising 
model [281430] . proteins [3T1435] . polymers [36], atomic clusters [37] and glasses [38l4i3] . 

In this work, we extend our previous results [22], which showed that quenched 
disorder lifts degeneracies in how the magnetic moments respond to a global driving 
field, allowing access to states that cannot be accessed in a perfect system. Here, we 
study how not only disorder but also driving field strength affects the accessibility of 
states and controllability of dynamics. 

The structure of this paper is as follows. In section [2] we describe our model 
system and outline the methods we use. In section [3] we discuss how the degree of 
a network node, that is, the number of links pointing into or out of it, relates to the 
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Figure 1. The geometry of the 4x4 array studied, (a) One of the two possible ground 
state tilings. The other ground state is obtained by flipping all spins, (b) One of the 
four possible polarised configurations, as would be obtained by applying a strong field 
in the +a; direction. 



energy of the configuration it represents. In section H] we study the number of states 
that can be reached from a polarised configuration and the reversibihty of dynamical 
transitions between those states. These results are built on in section [5l where we study 
the structure of the spin ice dynamics networks and discuss how small changes to the 
properties of the artificial spin ice can lead to large changes in dynamics. We give 
definitions throughout this paper of the network theoretic terminology and concepts 
used, but readers wishing for a more thorough introduction to network theory should 
see, for example, [Hl - HT] . 

2. Methods 

2.1. Energetics, dynamics and disorder 

In this work, we study a 4x4 spin ice array, with geometry shown in figurelH The 16 Ising 
spins of the array can take a total of 2^® = 65, 536 configurations. We show elsewhere [22] 
that this small, feasible to analyse, system has dynamics that are sufficiently similar to 
those of larger arrays that the results presented here are broadly relevant, even if small 
arrays are limited in their support of configurational features such as domain walls that 
are seen in larger systems. Larger arrays are inaccessible to our network analysis due to 
exponential growth in the number of configurational states (network nodes) with system 
size. For example, increasing the array size to even 5x5 spins increases the number 
of configurations - or equivalently, network nodes - to around 33, 000, 000. Of course, 
numerical simulations of dynamics of larger systems are effectively a sampling of their 
phase space networks, but we leave a rigorous interpretation of simulated dynamics in 
terms of networks for future work. 

We first describe our model for the energetics and dynamics of this system, and 
then describe how we construct a network representation of those dynamics. The island 
magnetic moments are Ising point dipoles, interacting so that the dipolar energy of 
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We set both the island magnetic moment M and the nearest-neighbour distance to 
unity, so that the nearest-neighbour couphng has strength 1.5, in units of /io/(4vr). This 
gives the energy scale of the island-island interactions. The total dipolar energy is given 
by (1/2) -^dip- ^ ground state spin configuration that minimises the total dipolar 
energy is shown in figure [TJ^a); the other ground state configuration is obtained by a 
global spin flip. 

The other field acting on each island is the external field h, which gives a Zeeman 
contribution to the total energy of spin i 



When a field strong enough to overcome the dipolar interactions is applied at 
approximately 45° to the island axes, the system's favoured state is a polarised 
configuration, in which all spins have the same projection onto the field. One of the 
four possible polarised configurations is shown in figure [TJ^b). Polarised configurations 
are of particular interest because they are experimentally reproducible. 

In addition to the island interactions, a second energy scale in the system is a 
barrier to island switching. We model this using the switching criterion 



component of the total field antiparallel to an island's magnetisation be greater than the 
island's intrinsic switching field hc^. A similar threshold-based model for switching has 
been used by other authors [7t[8|[TH[T^ . (We have studied other switching criteria, such 
as Stoner-Wohlfarth switching, in numerical simulations and find qualitatively similar 
dynamics.) In experimental systems, the island switching fields are usually designed 
to be larger than the dipolar fields so that configuration states can only change under 
external fields p. Here we set the mean switching field to he = 11.25, a value outside 
the range of dipolar coupling strengths. 

The dynamics of the system under fixed external field h consists of a series of single 
spin flips, determined by criterion (jl]). The set of all spins that satisfy (jlj) is calculated, 
then one is chosen uniformly at random and flipped. The set of spins satisfying (jl]) is re- 
calculated for the new conflguration, and again one is flipped at random. This process 
continues until no further spin flips are possible, at which point a 'flnal', stationary 
conflguration has been attained. Because spin flips are selected at random and the set 
of spins that can flip is recalculated at each step, more than one series of spin flips may 
be possible under application of the same fleld. In previous simulation studies [17], this 
meant that different simulation runs did not necessarily have the same outcome, and 
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we instead averaged over outcomes. In our network studies, we enumerate over all these 
possibilities, as described in section 12.21 

As seen in equation (jlj), the response of an island to the external driving field is 
controlled both by its interactions with other islands, and its intrinsic switching field, 
unless the external field is strong enough to overcome these. Inter-island interactions 
bias the response of the system towards low-energy states. At the same time, disorder, 
in the form of a distribution of island properties, introduces randomness in the response 
to fields. In this work, we focus exclusively on switching field disorder. We show 
elsewhere [21] that other types of disorder, such as disorder in island positions or 
orientations, have a similar effect on dynamics as switching field disorder does, and 
that all disorder can be characterised in terms of an effective switching field disorder. 
Based on that result, in this work we consider only switching field disorder. A direction 
for future study might be a study of disorder with correlations. 

We characterise the strength of disorder by a, the standard deviation of the 
switching fields for all islands. The switching fields are drawn from a nominally Gaussian 
distribution. However, because in each disorder realisation only 16 hc^ values are 
required, the mean of the generated pseudo-random numbers can deviate substantially 
from the nominal mean of 11.25. This change in mean switching field can have a 
significant effect on dynamics, so we 'correct' each hc^ value by —A, where A is the 
difference between the mean of the generated values and 11.25. Thus, for each disorder 
realisation studied, the mean switching field is exactly 11.25, the same as the switching 
field used in the absence of disorder, allowing meaningful comparison between disorder 
realisations. 

In our studies on field strength, we analyse three networks at each field strength: 
one for the perfect system, and two realisations of disorder. Both realisations are in the 
strong disorder regime of [21], with a = 2.05 and cr = 2.22 both larger than the scale of 
dipolar interactions. We focus on two similar disorder strengths in order to verify that 
networks representing systems with similar disorder strength have similar properties, 
regardless of field strength. 

2.2. Network construction 

The 4x4 Ising spin system we study has 2^^ = 65, 536 possible microscopic 
configurations. As mentioned in the Introduction, each configuration is a network node, 
and the number of nodes is fixed at 2^^ for all networks we study. For a given external 
field h, a link exists from node i to node / if the configuration corresponding to i can 
evolve into the configuration corresponding to / under h, via a cascade of spin fiips 
according to the dynamics described above. As we describe below, the network can be 
represented by a 2^^ x 2^^ matrix whose rows and columns represent nodes, and whose 
entries represent links (an 'adjacency matrix', in the language of network theory). 

The links of our networks are directed, that is, the existence of a link i — )■ / 
does not necessarily imply the existence of a link / — )■ i. The reason for this is that 
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dynamical transitions are not, in general, reversible: the athermal field- driven dynamics 
must involve transitions that lower the sum of dipolar and Zeeman energies, due to 
the nature of the energy barriers in the system. Directed networks appear in many 
other contexts, such as the world wide web [H], networks of corporate ownership and 
control [19], and the network representation of basins of mutually-reachable states in the 
random field Ising model |28] - [30] . In contrast, the networks used to describe, e.g, the 
six-vertex model [251 - I27] are undirected, because dynamics in those systems are taken 
to be reversible at the microscopic level. 

The set links of the network depend on the fields used to construct the network. 
For simplicity, in this work, each network we consider is constructed for a single field 
amplitude h and for field angles 6 = 0, 7r/128, 27r/128, . . .. This choice of field angles 
gives a network with properties approaching the expected finite limit for a continuously 
varying field angle, while remaining computationally tractable |22]. The network also 
depends on the disorder realisation used, with different realisations giving different 
networks. 

In order to enumerate all network links, we first determine all allowed single spin 
flips for each of the 2^^ configurations, for all fields {h,{6}). These can be stored as 
transition matrices T{h,6), where Tij{h,6) = 1 if configuration i can be transformed 
into configuration j by a single spin flip that is allowed under a field {h, 6) according to 
criterion (jlj), and Tij{h,6) = otherwise. Note that T{h,6) is not symmetric, because 
only flips that lower the system's total (dipolar and Zeeman) energy are allowed. 

The non-zero entries of give permitted transitions involving two spin flips. 
Similarly, three spin flips are described by T^, and so on. The non-zero entries of 
T* are not all equal to 1, so we apply the sign function (sign (a;) = 1 if x > and 

otherwise) to each element of T* to 'normalise' the matrix. Because there are 16 
spins in the system and spins can only flip once in a dynamical cascade, the maximum 
length of a sequence of spin flips is 16, and we must have sign(T*) = sign(T^^) for all 

1 > 16. We denote sign(T^^) as A{h,6), and Aij{h,6) = 1 if configuration i can evolve 
by a cascade of spin flips into configuration j, under the field {h,9), and Aij{h,9) = 
otherwise. We define the network adjacency matrix A{h) by A{h) = sign(^g A(/i, 6')), 
where the sum is over all 6 between and 2n. A describes the network representation 
of all possible dynamics under fields {h, {9}), and it is the properties of A for various 
field strengths and disorder realisations that will concern us in the rest of this paper. 
We emphasise that this method of network construction is an exact enumeration over 
all possible transitions between configuration states allowed for a given set of external 
fields. 

We also emphasise that each network link is active only for a range of field angles 
^min < < 6'max- Accordingly, the existence of a network path from a node i to a node 
j indicates that there exists one or more sequences of field angles that can be applied 
to a system prepared in state i to drive it into state j. A particular sequence of field 
angles, also known as a 'field protocol' in the literature, corresponds to the paths on 
the network that are generated by following links that are active for each angle in the 
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Figure 2. The (a) out degree and (b) in degree of 1000 randomly selected network 
nodes (spin configurations), plotted against their dipolar energy. Blue (dark grey) 
circles represent the perfect system, orange (light grey) squares represent the disordered 
system with a — 2.05. 



sequence in turn. Conversely, the network contains information about all possible field 
protocols. 

3. Node degrees and energetics 

In this section we discuss how an entirely local property of network nodes, namely their 
degree, can be related to spin ice physics. We also discuss briefly what the distribution 
of node degrees shows about the global topology of the network; however, we will see in 
subsequent sections that the degree distributions are insufficient to completely describe 
the network topology. 

The in- (out-) degree of a network node is the number of links pointing into (out 
of) it. In an undirected network, the two quantities are the same, but in a directed 
network they are different. The degree distribution 7\/'(/cin(out)) is the distribution of the 
number of nodes with degree fcin(out) • The in degree of node v is given by 

/Cjn = ^ ^ — v4^,^,, (5) 
i 

and the out degree is given by 

^out = ^ ^ Ayj A,^J,^, . (6) 

i 

Figure [2] shows the in- and out-degrees of 1000 randomly selected nodes, vs the 
energy of the spin configurations they represent, for networks describing a perfect and 
a disordered system (with a = 2.05) subject to a field of amplitude h = 11.5. For both 
perfect and disordered systems, high-energy configurations correspond to nodes with 
low in-degree and high out-degree, and low-energy configurations correspond to nodes 
with a higher in-degree and lower out-degree. 

Physically, a node with high out-degree represents a configuration that can 'decay' 
into many others when fields are applied, while a node with low out-degree represents a 
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'stable' configuration. 'Stability' here refers to how a configuration can be modified by 
an external field, as measured by, for example, the probability that a field with direction 
chosen uniformly at random is able to drive the system to a new configuration. Although 
the stability of a configuration is not completely correlated with its energy and the out- 
degree vs energy data displays some spread, there is a clear relationship between energy 
and out-degree, as seen in figure [2](a). This can be understood by considering the 
extremes of high and low energy: At the field strength studied here {h = 11.5), all 
spins in a configuration that maximises dipolar energy can be flipped by an external 
field even when moderate disorder is present, so there are many pathways out of that 
configuration. On the other hand, the ground state of the system is stable even under 
moderate disorder. 

The energy dependence of the in-degree has a wider spread, especially when in- 
degree is low. This is because the barriers to 'entering' a state are topological as well 
as energetic. A low in-degree may correspond to a state with high energy or it may 
correspond to a low energy state with antiferromagnetic ordering that is 'hard' to access 
with a global external field, which tends to create ferromagnetic ordering. We discuss 
these ideas in relation to the ground state of the system elsewhere [1]. Nevertheless, it is 
the case that only low-energy configurations have high in-degree, as seen in figure [2]^b). 

Node degrees can also give information about the global network topology, via the 
degree distribution jV{k). For example, many real- world networks such as the world wide 
web and networks of scientific citation display power-law degree distributions P^IHS] . On 
the other hand, in Erdos-Renyi random graphs [50], where an undirected link between 
any pair of nodes is present with probability p and absent with probability I — p, 
the degree distribution is Poissonian. Networks describing the phase space of other 
frustrated spin models and lattice gas models have been shown to have Gaussian degree 
distributions [25H27]. 

Directed networks can be described by three distributions: the joint in- and out- 
degree distribution, which gives the probability that a randomly selected node has 
in-degree fcin and out-degree fcout, and the two separate degree distributions, which 
are obtained by integrating the joint distribution. In figure |3l we plot these three 
distributions for the network describing an undisordered spin ice at /i = 11. In contrast 
to other frustrated systems [251 - I27] . the spin ice networks do not have a Gaussian 
degree distribution. Instead, the distributions are clearly asymmetric, and the separate 
distributions for in- and out-degree are different. As expected from our discussion above 
about the relationship between configuration energy and node degree, there is a tendency 
for low in-degree nodes to have a higher out-degree, and vice-versa. However, we will 
see in section [5] that the degree distribution is inadequate to completely characterise the 
network. 
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Figure 3. Degree distributions for a perfect system, subject to a field ft, = 11. (a) 
The joint in- and out-degree distribution, for fcin < 50 and fcout < 30. (b) The separate 
out-degree distribution, (c) The separate in-degree distribution. Disordered systems 
have similar degree distributions. 



4. Accessibility of states and reversibility of dynamics 

Many potential technological applications of artificial spin ices and related systems 
depend on the answers to two questions. How many distinct configurational states are 
available, and can those states be accessed reliably? These questions are particularly 
important for applications relating to information storage: although a system of Ising 
spins has in principle 2^ configurations and can store N bits of information, the effective 
information capacity is lower if only some fraction of those configurations can be realised. 
In this section, we see that network tools are ideal for studying these questions, and for 
uncovering effects of disorder and field strength on the accessibility of states. 

As already mentioned in section I2TT1 the only magnetisation configurations currently 
known to be exactly experimentally reproducible are the four polarised configurations 
(see figure mb)). Because these configurations can be reliably obtained, in almost all 
experimental and simulation studies to date, the initial state of the system has been 
polarised [Il|5l|10llI3[Tl[l5l|l71ll9l|20l[5T]. It is therefore of particular interest to study 
the accessibility of states from these configurations. 

We focus on states that can be accessed using a sequence of fields with fixed 
amplitude. As mentioned in section 12.21 such a sequence of fields corresponds to a 
'walk' on the network, in which at each step only a link that is active at that field angle 
is followed. There may be more than one walk corresponding to a sequence of fields, 
since there may be more than one link active for a given field angle. A random field 
protocol, in which the field angle is selected uniformly at random from [0, 2tt) at each 
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Figure 4. (a) The fraction of configurations that can be reached from an initial 
polarised configuration, vs applied field amplitude, (b, c) The fraction of configurations 
that can be reached from an initial polarised configuration, vs disorder strength 
(standard deviation in island switching fields), for an applied field of strength h = 11.5 
(b) and h ~ 11.75 (c). For h = 11.5, there is a jump in the number of accessible states 
when disorder is turned on, but for h = 11.75, the number of accessible states increases 
smoothly. Each point represents a single disorder realisation. In all figures, the initial 
configuration is counted in the number of accessible states. 



step of the sequence, is essentially a random walk on the network 

The total set of nodes that can be reached from an initial node v, following network 
paths of any length, is given by the fixed point of repeatedly multiplying the unit vector 
V (all entries zero, except for non-zero entry v) by the adjacency matrix A. Each state 
can be accessed via one or more field protocols, but a given field protocol may not be 
able to access all of them. We return to this question of ergodicity below. 

Figure Hl^a) shows how the fraction of configurations that can be reached from the 
+x polarised configuration depends on field strength, for both the perfect system and 
two disorder realisations (cr = 2.05 and a = 2.22). For weak fields, the +x polarised 
configuration is stable against applied field and no other states are accessible from it. 
However, for optimal field values, approximately 10% of all configurations can be reached 
from the +x polarised configuration. In the very high field limit (not shown), applied 
fields always polarise the system, so from the +x polarised configuration there are 3 
other states accessible, namely, the other polarised states. 

Although the curves for the perfect and disordered systems in figure IH^a) have 
similar form, for applied field strengths near the mean island switching field [h^ = 11.25) 
the difference between perfect and disordered systems is large. This is in agreement with 
previous results [22]. Figure Ht^b) shows how the number of states accessible from the 
+x polarised configuration depends on disorder strength, for h = 11.5. There is a 
jump in the number of accessible states when disorder is turned on. This is because 
two of the configurations that can be reached from the +x polarised configuration have 
spins that, in the perfect system, require an external field of 11.74 to switch. A small 
disorder-induced decrease in the switching barrier for these spins allows them to flip at 
h = 11.5, opening new dynamical pathways. This interpretation is confirmed by the 

I The equivalence between the simulation of a random field protocol and a random walk on the network 
is not exact, because the network is unweighted and the probability a link is followed in a random walk 
depends only on the number of links out of a node, not the range of field angles over which they are 
active. 



A network model for field and quenched disorder effects in artificial spin ice. 



11 




Figure 5. (a) An example directed network containing 4 nodes. The nodes A, B, 
C form a strongly connected component, but D is not a member of the component 
because there is no path from D to the other nodes, (b) An example of how a node 
may be accessible, but not reproducibly so. The highlighted node can be reached by 
applying a field of angle 9 to node A, or a field of angle a to node B, but both of these 
field angles also activate links to other nodes, so it is not possible to construct a field 
sequence that is guaranteed to pass through the highlighted node. 

network for h = 11.75. For that field strength, the fraction of states accessible from 
the +x polarised configuration in the undisordered system is ~ 10~^, and increases 
smoothly with disorder strength, as shown in Figure |4](c). The large impact of small 
disorder-induced changes to the system is a recurring theme in this work. 

The number of states that can be reached from an initial polarised state is a starting 
point for describing dynamics, but this quantity does not give a complete picture. One 
question it gives little information about is that of ergodicity - once a transition has 
been made from the polarised state to another state, what further transitions can be 
made? Is it still possible to access all of the other configurations that are accessible 
from the polarised configuration? Can transitions be reversed? This is important for 
information storage applications, where it is important that the state of the system can 
be easily and reliably 're-written'. 

A useful network theoretic concept here is that of strongly connected components 
(SCCs). An sec of a directed network is a set of nodes for which paths exist between 
every pair of nodes, taking the directions of the links into account. For example, in the 
network shown in figure [5](a), the nodes A, B, C form an SCC, but D is not a member 
of the component because there is no path from D to the other nodes. We determine the 
SCCs of a network using the algorithm in ^52j, as implemented by the software package 
Mathematica. 

Dynamics within an SCC are reversible, provided the correct field protocol is 
applied. In terms of networks, it is possible to travel along network links from any 
node of the SCC to any other node of the SCC; in terms of artificial spin ice dynamics 
this means that for any configuration in the SCC, there exists a sequence of fields to 
drive the system from that configuration to any other configuration in the SCC, and 
back again. 

The sizes of two SCCs in particular are of interest: the SCC that contains the +x 
polarised configuration (or one of the other three polarised configurations: this choice is 
arbitrary), and the size of the largest SCC in the network. Figure [6](a) shows these sizes, 
plotted against field strength, for the system without disorder and the two disordered 
systems. For the system without disorder, the +x polarised configuration is always in 
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Figure 6. (a) Fraction of all nodes that are in the largest strongly connected 
component (SCC) and the SCC that contains the +x polarised configuration, as 
a function of field strength. For the system without disorder, the +x polarised 
configuration is always in the largest SCC. (b)Fraction of all nodes that are in the 
largest SCC, vs disorder strength (standard deviation in island switching fields). The 
applied field strength is h = 11.5. Each point represents a single disorder realisation. 

the largest SCC. 

Below a threshold field, all SCCs have size 1, that is, no dynamical transitions can 
be reversed. Above the threshold, the size of the largest SCC grows by three orders of 
magnitude to take in approximately ten percent of all nodes, before decreasing again 
for large field strengths. In the limit of very strong fields the largest SCC has size 4, 
and consists of the four polarised configurations. This limit holds for both perfect and 
disordered systems because strong external fields overcome disorder. When disorder is 
present, the +x polarised configuration is not always in the largest SCC, but it is when 
the external field is sufficiently strong. 

Comparison of figures H] and [6] indicates the strong correlation between the number 
of accessible states and SCC size. In general, the number of states accessible from a 
node must be greater than or equal to the size of its SCC. Figure[5](a) illustrates a simple 
example of this: the number of nodes reachable from A is 4, while it is in an SCC of 
size 3. In fact, in the perfect system, for h > 12, all four polarised configurations are 
in the largest SCC and for h > 14 all states that can be reached from the +x polarised 
configuration are in the same SCC. Similar results hold for the disordered systems. 

The existence of large SCCs that include the polarised configurations implies that 
for correctly-tuned fields, the information storage capacity of the artificial spin ice is 
maximised, with several thousand configurations accessible from one another, making it 
possible to 'write' a configuration and then 'rewrite' a new configuration by applying a 
suitable sequence of fields. However, one should be careful for two reasons. First, if the 
number of states accessible from a given starting state (e.g, a polarised configuration) is 
larger than the SCC size, there will be dynamical 'dead ends', that is, states that can be 
entered but not exited. Second, the existence of a path into a node does not guarantee 
that it is possible to reliably access that configuration. This is illustrated in figure [5t^b), 
where the highlighted node can be accessed, but it is not possible to construct a field 
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sequence that is guaranteed to pass through it. A more detailed study of these two 
points is a topic for future work. 

We close this section by commenting on the effect of disorder. While the general 
trends for SCC sizes as a function of field strength are the same for disordered and 
perfect systems, for applied fields close to the mean island switching field of 11.25, the 
difference in SCC size between perfect and disordered systems is substantial: around 
two orders of magnitude for h = 11. Figure |6]^b) shows the size of the largest SCC for 
a range of disorder realisations, at an applied field strength of h = 11.5. The size of 
the largest SCC increases continuously with disorder strength, although, as might be 
expected, the spread in values for different realisations of strong disorder is substantial. 

Although the size of the largest SCC grows with disorder strength in the regime we 
study, we note that for a given array size and field strength, there is an upper bound on 
the number of accessible states, which is proportional to 2^', where A^' is the number 
of spins that are not pinned and that do not always align with the external field, that 
is, the number of spins that are not 'frozen' by disorder. This quantity decreases with 
disorder strength. This argument suggests there should be an optimal disorder strength 
for accessibility of states, but the optimal value depends on array size and field strength, 
rather than being universal to artificial spin ices. 

5. Network structure, 'rewiring', and control 

We saw in section[3]that the energetics of artificial spin ice plays a key role in determining 
network properties on a local, i.e: node, level. Similarly, in section HI we have seen 
how the global properties of the network relate to dynamics, which can be observed 
in simulation and are experimentally testable. In this section, we comment on the 
relationship between the local and global scales of the network, showing that purely 
local network properties are insufficient to predict the global structure and that instead, 
correlations exist in the network that are determined by the dynamics of artificial spin 
ice. We then show how the insight this gives about the network structure can be used 
to explore the possibility of controlling artificial spin ices. 

We first demonstrate the existence of correlations in the artificial spin ice dynamics 
networks. We do this by comparing the spin ice networks with two other type of 
networks that retain some properties of the spin ice dynamics networks but are otherwise 
randomised. Comparison with randomised networks is frequently used to reveal the 
structure of real- world networks, see, for example, [53H56] . 

The first type of randomised network consists of uncorrelated random networks |50] , 
in which any link i — )■ / exists with equal probability, p = n(links)/2^^ where n(links) 
is the number of links and 2^^ is the total number of possible links between 2^^ nodes. 
The uncorrelated random networks represent a dynamics in which the possibility to pass 
from any state to any other state does not depend on energy at all. In particular, the 
relationship between the energy of a configuration and the degree of its network node 
is lost. 
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Figure 7. The fraction of nodes in the largest strongly connected component (SCC), 
vs number of network links for spin ice dynamics networks and 'test' networks with 
fewer correlations. In the 'uncorrelated' test network, any link i — > / exists with 
equal probability. In the 'randomised' test networks, the joint in- and out-degree 
distributions are the same as those of the undisordered spin ice network, but the links 
themselves are random. 



A more sophisticated approach to constructing random networks is to preserve this 
relationship, by creating 'maximally random' networks consistent with a given joint in- 
and out-degree distribution [57l|58] (see section |3]). Under such a scheme, high-energy, 
unstable configurations have many links out, and low-energy, stable configurations have 
few links out, for example. This gives a network that, at least locally, represents a 
dynamics much more closely approximating the actual artificial spin ice dynamics. 

The comparison between these networks and the spin ice dynamics networks is 
illustrated in figure [TJ where we plot the fraction of nodes in the largest SCC (strongly 
connected component; see section H]) against the number of links in the networl|§|. For 
the spin ice dynamics networks, the size of the largest SCC is not a single- valued function 
of the number of links. However, as seen in figure [6]^a), it is a single- valued function of 
field strength, which parameterises the curve, and the 'doubling back' of the SCC size 
vs the number of links occurs because the number of links has a peak near h = 17. 

The randomised networks have very different global properties to the spin ice 
dynamics networks. This is because the connections in the spin ice dynamics networks 
are not simply dictated by the degree distribution, which in turn is because the dynamics 
are not simply dictated by the energies of states but also by the barriers between them. 
While all networks show an increase in the size of the largest SCC with number of 
network links, the spin ice djTiamics networks have a much reduced tendency towards 
large SCCs than the randomised networks do. SCC sizes are commonly used as a 
measure of percolation in directed networks [57 1 159] - [63] . but, unlike the randomised 
networks, the spin ice dynamics networks never fully percolate. This may indicate a 

§ A note on definitions: If nodes A and B are linked in both directions, that is, {A — > B) and {B ^ A), 
we count each directed link separately, to give two links. Unless otherwise specified, we do not count 
self-links, that is links pointing from a node to itself. Under this definition, the number of links is given 

by J2^^J^^J. 
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high level of clustering in the network |64j . 

An intriguing feature of figure [7] is the dramatic increase in the size of the largest 
sec in the artificial spin ice networks caused by a relatively small increase in the 
number of network links. When the SCC is enlarged by tuning the field strength, the 
corresponding increase in the number of links is approximately ten-fold: from 622, 896 
to 6, 180, 266 in the field window over which the largest SCC size in the perfect system 
increases from 1 node (10~^ of all nodes) to 8% of all nodes. When disorder is used 
to increase the SCC size, the change in the number of links is even smaller: at fixed 
field strength h = 11.5, disorder increases the number of links from 736, 720 in a perfect 
system to 1, 060, 814 when the disorder strength is a = 3.3 - a change that is associated 
with the largest strongly connected component size growing from 3 nodes (10~^ of all 
nodes) to 19% of all nodes. 

The fact that a relatively small change in the links of the network can alter 
the connectedness of nodes so dramatically suggests that for fields close to the mean 
switching field of 11.25 the perfect system is 'almost' well connected. This notion is 
supported by the jump in the number of states accessible from a polarised configuration 
at /i = 11.5 when disorder is turned on, shown in figure Hl^b), which we have already seen 
is caused by disorder allowing spins to flip that in the perfect system have a switching 
barrier slightly higher than the applied field. 

The possibility for small changes to the system to have large effects on the 
accessibility of states finds application in the control of artificial spin ice. For example, 
in experimental studies of field-driven reversal in artificial kagome spin ice, islands were 
deliberately modified in order to serve as 'start' and 'stop' sites for avalanches of spin 
flips [8j, and simulations of a colloidal model for artificial spin ice reveal that using 
different barrier heights for different sublattices leads to a rich array of stable states 
that are different to those seen when all barriers are the same [65]. As an alternative, 
one might imagine a system where a small number of macrospins are controlled directly 
via, for example, current-driving switching. In a network picture, such modifications of 
the spin ice system are equivalent to deliberately creating and removing certain links. 
As seen already in this section, such re- wiring can have a dramatic effect. 

This effect is illustrated by the networks shown in figure [HI Network (a) is the 
network of states accessible from the +x polarised configuration (the large red node) 
for a perfect system, at a field strength of h = 11.5. There are two other configurations 
accessible via a single field application, and a further two accessible if a second field 
is applied, making a total of 5 nodes. In contrast, network (b) is the network of 
configurations accessible from the same initial configuration, for the same field strength, 
but in a system where the lower-left corner island can be flipped independently of the 
others. This network contains 128 nodes. In other words, the ability to separately 
control a single spin yields an order-of-magnitude increase in the number of states 
accessible from the +x polarised configuration. 

These preliminary results demonstrate the value of the network picture of artificial 
spin ice dynamics for studying these problems. Future work might take advantage 
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Figure 8. The network of states accessible from tire +x polarised configuration (the 
large red node), at a field strength oi h = 11.5, for (a) a perfect system, and (b) a 
system where the lower-left corner island can be flipped independently of the others. 

of other network properties. For example, in studies of how epidemics spread on 
networks, tools have been developed to determine which nodes are most important 
in determining the properties of transport on the network Since field- induced 

dynamics are essentially the same as network transport in this picture, applying these 
tools to networks describing spin ices may offer a way to determine how to modify the 
spin ice to control dynamics as desired. 

6. Conclusion and outlook 

To summarise: we have shown that a network model for the dynamics of artificial spin 
ice provides a means of quantifying how applied field strength and quenched disorder 
affect the system's behaviour. Increasing disorder strength and tuning the applied field 
increase the number of states accessible to field-driven dynamics and the reversibility 
of dynamical transitions, both of which are important for potential applications. The 
changes in dynamics are caused by a 're-wiring' of the network that involves relatively 
few links, suggesting that the highly restricted dynamics of a perfect system subject to 




A network model for field and quenched disorder effects in artificial spin ice. 17 

a sub-optimal field are caused by a small number of dynamical pathways being blocked. 
We have shown that, indeed, a small change to the artificial spin ice system unblocks 
these pathways and allows many new configurations to be accessed. We have also shown 
that the degree of a network node, a local property of the network, can also be related 
to the physics of the system, via the correlation between node degree and the energy of 
the configuration it represents. 

It would be interesting to apply the network tools we have developed here to 
understand other problems in artificial spin ice. For example, the square ice studied here 
is only partially frustrated, and has a well-defined ground state even in the limit of short- 
range dipolar interactions. On the other hand, artificial kagome ice pl[3] or a proposed 
square ice with sublattice height offsets |TT1|66] are fully frustrated with extensive 
degeneracy, at least when interactions are short-ranged [TTll67[l68j . Experimental studies 
have also been made comparing the demagnetisation of frustrated and unfrustrated 
small clusters of islands [69] . Comparing the network properties of a fully frustrated ice 
with the partially frustrated square ice may give clues to the role of frustration in the 
dynamics. 

There are also open avenues of enquiry related to the network structures themselves. 
For example, in studies of frustrated triangular antiferromagnets and the six vertex 
model, Peng et al. [27j find fractal structures in the phase space networks which they 
suggest may be a signature of 'long-range interactions, correlations, or boundary effects 
in real space'. The artificial spin ice systems studied here exhibit long range interactions, 
but whether their phase space networks have fractal structure is not yet known. 
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